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Abstract. We investigate partition- function zeros of the many-body interacting 
spherical spin glass, the so-called p-spin spherical model, with respect to the complex 
temperature in the thermodynamic limit. We use the replica method and extend the 
procedure of the replica symmetry breaking ansatz to be applicable in the complex- 
parameter case. We derive the phase diagrams in the complex-temperature plane and 
calculate the density of zeros in each phase. Near the imaginary axis away from the 
origin, there is a replica symmetric phase having a large density. On the other hand, 
we observe no density in the spin-glass phases, irrespective of the replica symmetry 
breaking. We speculate that this suggests the absence of the temperature chaos. To 
confirm this, we investigate the multiple many-body interacting case which is known 
to exhibit the chaos effect. The result shows that the density of zeros actually takes 
finite values in the spin-glass phase, even on the real axis. These observations indicate 
that the density of zeros is more closely connected to the chaos effect than the replica 
symmetry breaking. 
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1. Introduction 

Phase transitions and critical phenomena have been a central problem in statistical 
physics for decades. After many pioneering works, it was revealed that phase transitions 
can be identified as singularities of the free energy, and several approaches to capture 
them have also been investigated for a long time [U, Ej |3l HJ [5] . The theory of partition- 
function zeros invented by Yang and Lee [3l H] is one of such approaches and offers 
a novel, and simple, picture of phase transitions. They proved that the free energy is 
analytic in a region where there are no zeros of the partition function, and hence there is 
no phase transition in that region. Besides, using Ising ferromagnets, they demonstrated 
that the phase transitions of the models become clearly visible by the zeros. Their work 
was followed by many other researchers and was applied to various situations [5l |6] 

Spin glass (SG) is known to show nontrivial phase transitions and critical 
phenomena and has been studied for a long time O |8]. According to the standard 
description of SGs, a SG system acquires a multi-valley structure in the free energy 
landscape at low temperatures. Some peculiar properties associated with SG transitions, 
such as strong hysteresis and rejuvenation-memory effect, are explained on the basis of 
this rugged landscape. In the mean-field level, each valley of the free energy is separated 
by infinitely-high free-energy barriers and is called a pure state. Each pure state 
corresponds to a thermodynamic phase. This provides a speculation that a sequence of 
phase transitions in a sense can occur in SG phases where the dominant part of pure 
states can vary as we change external parameters such as temperature. Unfortunately, 
it is difficult to directly examine this kind of transitions in the mean-field solution given 
by Parisi [9l [10] . This motivates us to use another approach, the zeros theory by Yang 
and Lee. 

Zeros of SGs have also been investigated for a fairly long time [HI [121 [131 [IH [13 
[T6l [T71 [T8l [To] , but until very recently the reasonable solution in the thermodynamic 
limit was not obtained even for mean-field models [201 [21], except for a simple model 
named the random energy model (REM) [161 HH]- Observing these few solutions 
in the thermodynamic limit [161 [201 [21], we find a tendency that the distributions 
of zeros are closely related to the step number of the replica symmetry breaking 
(RSB) in the Parisi scheme. For Bethe SGs exhibiting the full-step RSB (FRSB), 
the zeros tend to densely distribute around the real axes of the temperature and 
uniform field below the critical points [20]. This suggests that a certain type of phase 
transitions occur everywhere in the FRSB phase, which supports the above speculation. 
This type of transitions can be possibly interpreted as the temperature/field chaos 
meaning that the spin configuration drastically changes as the temperature/field slightly 
varies [22l [231 [22 [251 [261 [23 |28l [291 [30]. On the other hand, for a family of REMs 
exhibiting the one-step RSB (IRSB), no zeros in the complex temperature plane, 
or very few zeros in the complex field plane, exist in the internal region of the SG 
phases [161 [HI [21]. These contrasting results possibly reflect the difference between 
the IRSB and FRSB, or may be just due to the peculiarity of the REMs being over- 
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simplified models. To make this point clear, we need to analyze distributions of zeros 
in more realistic SG models, which is the main purpose of this paper. 

Although the validity of the RSB picture is questioned in finite-dimensional 
SGs[3ll|32l|33l[Ml[35l[Ml[371|38l|39lll^ temperature chaos is considered 

to exist in those systems [321 El HHl SSI SH SHI US]. Hence, it will be helpful to 
reveal the relation between the chaos effect and the density of zeros (DOZ). Up to 
the present, we have no fully-reliable result about the zeros of finite-dimensional SGs 
in the thermodynamic limit, though they are intensively studied [HI [I3l [HI [I5]. We 
expect that some clear knowledge about the DOZ of SGs, even in the mean-field level, 
can be a help to improve this situation. 

In this paper, we investigate the distribution of zeros of the many-body interacting 
spherical SGs [50l [51]. These are more natural than the REMs in that the phase-space 
decomposition into many pure states occurs as the temperature changes. Besides, it 
shows both the replica symmetric (RS) SG phase and the IRSB-SG phase depending 
on a parameter p being the number of interacting spins. Moreover, it is easy to control 
the temperature chaos of this model. In the case of the single p-body interaction, there 
is no temperature chaos. However, in the multiple {p + r)-body interacting case it is 
known that the temperature chaos occurs [30]. These properties are quite useful to 
investigate the relations among the DOZ, the RSB and the chaos effect. 

To derive the DOZ, we use the formulation invented in [2r|. In that formulation, 
we employ the replica method and generalize the Parisi scheme to be applicable in 
the complex-parameter cases. Three different types of overlaps between replicas are 
introduced. The physical interpretation of the overlaps is also one of the results in this 
paper. Although we concentrate only on the zeros in the complex temperature plane, 
our formulation can be applied to the complex field or other physical parameters. 

This paper is organized as follows. In the next section, we start from a brief 
introduction of partition-function zeros. The replica-based formulation to assess DOZs 
is also explained in this section. In section [3l we introduce the spherical SG model and 
derive the saddle-point equations to calculate the DOZ in the replica formulation. The 
RS and IRSB solutions are investigated and the physical significances are discussed. 
In section [H we present the phase diagrams in the complex temperature plane and the 
values of the DOZ in each phase. Last section is devoted to conclusion. 

2. Formulation 

2.1. Partition- function zeros 

Since a partition function of a finite size system is generally analytic with respect to a 
physical parameter y, we can reasonably assume that the partition function of the size 
can be factorized into a product form 





(1) 



j 
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where {y^-'-'} are zeros of the partition function and generally complex y^-'-' = y[^^ + iy2\ 
where i is the imaginary unit. We assume that C is an analytic and irrelevant factor. 
Hence, the free energy density of the system / = —{Nl3)~^\n Z is written by 

- = C + 5^ ^ \n{y - = C + J dz,dz2p{z,, z,) \n{y - z), (2) 

where z = zi + iz2 and we define the DOZ p{zi, Z2) as 

p{zuZ2) = ^J2^{z-y^^^ (3) 
j 

Since C is analytic, singularities of the free energy are characterized by the DOZ only, 
which motivates us to investigate the DOZ. 
The delta function can be rewritten as 

which is the same relation as the one between a point charge and an electrostatic 
potential in electrostatics in two dimension. This relation leads to 

This formula becomes a base of the following discussion. 

Note that (I5j) is the Poisson equation, which means that the relation between 
the generating function g{y) and the DOZ p{y) is compared to the one between the 
electrostatic potential and the charge density in two dimension. The one-dimensionally- 
distributed charge density is evaluated from the discontinuity of the electric field. 
Similarly, the one-dimensionally-distributed DOZ appearing on phase boundaries can 
be assessed by the difference of first derivatives of g{y) in the adjacent phases. Based on 
this analogy, we can derive the following formula for the one- dimensional DOZ pid{yi, 2/2) 
on a phase boundary represented by b{yi,y2) = 0: 

f \ ^ \(dgi dg2\ db (dgi dg2\ dh\ 

where gi{y) and g2{y) are the generating functions of the adjacent phases. The function 
h is defined such that h is positive (negative) in the phase 1 (2) and becomes zero on 
the boundary given by gi = g2- In most cases, it can be chosen as b = gi — g2. 



2.2. Zeros of random systems and the replica method 

For random systems such as SGs, the DOZ fluctuates from sample to sample. In the 
thermodynamic limit, we can expect that the typical DOZ converges to the averaged 
one. This requires to take a difficult average of the logarithm as [ln|Z|], where the 
brackets [(■ ■ ■)] denote the average over the quenched randomness. The replica method 
bypasses this problem by using the identity 

g{y) = h\n\Z{y)\] = lim -1- ln[|Z(y)p"] ^ \im^<P{y,n). (7) 
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Once we obtain (j){y,n), we can calculate the DOZ from 4>{y,n) through and ([7]). 
Unfortunately, it is still difficult to treat the nth power for real n. To avoid this difficulty, 
we first assume that the exponent n is an integer and evaluate in that condition. 

After that, we take the limit n — )■ by utilizing the analytic continuation from integer 
to real. 

This standard prescription of the replica method has some delicate problems in 
taking n — )■ limit. In some cases, a naive analytic continuation (RS solution) leads 
to an incorrect result, and the RSB solution is required. The RSB takes the rugged 
landscape of the free energy into account, which is essential to consider SG systems. 
In the present formulation, the RSB is implemented as an ansatz in the overlap matrix 
among n replicas as usual. In our formulation to calculate = [(ZZ*)"], we have 

three types of overlaps: the usual overlap {q} among n replicas of Z, {q'} among n 
replicas of Z*, and the inter-overlap {q} between n replicas of Z and those of Z*. 
Hence, we need some modifications in the RS and RSB ansatz to treat this extended 
overlap matrix. The detailed discussion about this point is presented in the next section 
after constructing the replica solution of the spherical SG. 

Before closing this section, we mention a physical consequence of the inter-overlap 
{q}. If {q} vanishes, the generating function decouples as g{y) = [InZ + \nZ*]/N and 
the DOZ inevitably vanishes. However, the reverse is not necessarily true. Even when 
the inter-overlap takes a finite value, the DOZ can become zero. An actual example is 
shown in section HJ We also mention that the replica method to calculate the zeros is 
very similar to that to find the chaos effect |26l ETJ EH] EHl 152] • In both calculations, the 
replica space is doubled to find the nontrivial effects of inter-overlaps. 

3. Replica analysis of zeros of the spherical spin glass 

The Hamiltonian of the p-body interacting spherical SG is given by 



where the spin Si takes continuous values under the spherical constraint ^ ■ Sf = N, 
and the interaction Ji^...ip is drawn from Gaussian with the variance J'^p\/2N^~^ 



In this paper, we also treat the multiple {p + r)-body interacting case, but below we 
explain our formulation on the single p-body interacting case. This is for the simplicity 
of the notation, but the analysis of the multiple-body case is essentially the same as the 
single one and hence the generalization is straightforward. 

To assess the zeros in the complex temperature plane, we calculate [|Zp"] by using 
the replica method as noted in the previous section. We can write [|Zp"] = [Z"'{Z*)"'] 




(8) 



n<---<ip 
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|2ni 



a=l 



ii<---<ip \ a=l 



n / N 



X 



(10) 



a=l \i=l 



J=l 



where Tr means the integration over all the spin variables. The spherical constraint is 
expressed in the delta functions. The average [(• ■ ■)] can be easily performed 



exp 



ii<---<ip 



a=l 



a=l 



exp N 



a,b 



a,b 



2 j2 
2 



Tab ( ' 



a,b 



where we put 



Qab = E^ S^S^, q^fj = — 5'j°'S'/, qab = — 'S'f 5*^ 

i=l i=l 4=1 

Note that subleading terms are omitted by using the following relation 



pi 



(11) 



(12) 



(13) 



Ef Qa Qb \ ( na nb ^ 

ii<-<ip 

Let us express the relation ( fT2l) by delta functions and introduce the overlaps 
{labyQaby'lab} as integration variables into (fTT]) . Besides, we also rewrite the delta 
functions by the Fourier expressions as 



(14) 



The factor A is arbitrary and chosen to make the following calculations simple. Hence, 
f lTU]) reads 



|2ni 



J]^ dqabdq'^fjdqabdXabdx'abdXab exp 



where 



E ~ ^ E (^«''^«'' + ^afe^ab + "^XabQab) + lu Tr 6^ l , 
a,b a,b 



L = IY1 {xabS'S' + XabS'^S'' + 2XabS''S'' 
a,b 



(15) 



X X 
XT X' 



s 

S' 



(16) 
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The matrix X has the element Xab = Xab, and X' and X are defined in the same 
way. The matrix X'^ represents the transposition of X. The spherical constraint is also 
expressed by the Fourier expression and is absorbed into the diagonal part of X and X'. 
The factor L is quadratic and the spin integration can be easily performed. The result 
is 

InTre"^ = ^Trln { -2n I ^-^^ } . (17) 




Note that Tr in the right-hand side, and henceforth, denotes the trace of the matrix. 
In the thermodynamic limit, we can use the saddle-point method to evaluate ( !T5|) . 
Substituting ( !T7|) into ( ITSl) and taking the saddle-point condition with respect to X, X' 
and X, we find 



X X \ ^ Q Q 

X^ X' [ gT g/ 



The matrix Q is defined as Qab = Qab, and Q' and Q are defined in the same way. 
Note that the real part of Wq should be positive for the convergence of the spin 
trace performed in ( fTTj) . Summarizing the above transformations, we get 0(/3,n) = 
ln[|Z|2"]/X as 

a,b 

+ ^TrlniyQ + n(l + ln27r). (19) 

For further calculations, we need some ansatz with respect to Wq. Note that the 
diagonal part of Wq, qaa and is equal to unity qaa = Qaa = ^- This is the consequence 
of the spherical constraint, which can be easily seen in ( !T2|) . We also have the condition 
that the matrix Wq is symmetric as Wq = Wq. 



3.1. RS ansatz 

We start from the simplest RS case. Due to the additional overlaps coming from the 
complex parameter, even the RS solution has a nontrivial form. We first illustrate this 
point. 

Let us refer to the analysis of the REM given in [53]. In that analysis, there are 
two types of RS solutions even in the real-parameter case. They are distinguished by 
the way of partitioning the n replicas into the spin states. In one way, all the replicas 
are distributed into different states, and in the other one all the replicas are in a single 
state. Fortunately, these solutions are parameterized by a single overlap q and we need 
not to distinguish them in the form of the overlap matrix. The first way gives the 
paramagnetic solution g = and the other yields the SG solution q = I. 

To construct the correct RS form of the overlap matrix in the complex-parameter 
case, we need to generalize the above replica-partitioning way. Let us extend the 
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interpretation of the "state" into the pure statqj. We assume that either of the above 
two-types solutions is reahzed if we focus only on the original n replicas Z", or only on 
its conjugate (Z*)". Considering the symmetry between and (Z*)", we obtain four 
types of RS solutions, which are schematically expressed in figure [1] Then, under this 
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Figure 1. Four types of RS solutions. Each ball represents a replica and each box a 
pure state. The upper half of the boxes corresponds to the phase space of the original 
partition function Z, and the lower one is the counterpart of Z*. We assume that the 
number of pure states is larger than n e N. 



ansatz, the overlap matrices can be parameterized as 
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In each matrix, all the off-diagonal elements are filled in by the same parameter. We 
note that any permutation of Q also becomes a solution. This is because the indices 
of the original replicas and those of the conjugate ones (Z*)" can be independently 
chosen. We here present the simplest form. This solution includes that of REM with 
complex parameters [2T] . 

Two parameters of Q, qi and go, are needed to express the P2 solution and are the 
consequence of introduction of the complex parameter. Actually, the P2 solution gives 
two-dimensional distribution of zeros, which is a property distinguished from the other 
RS phases being characterized by qi = %. 

Using this RS ansatz, we next calculate n). We can write the first term in ( |T9l) 

Y.€, = n + n{n-l)q^, (21) 

a,b 

X This interpretation includes the REM case, since a pure state of the REM coincides with a spin state. 
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and the second and third ones as well. To calculate the Tr term in f ll9p . we need to 
diagonalize Wq. The calculations are simple and omitted here. The result is 

TrlniyQ = lnAi + (n-l)lnAo, (22) 

where 

Ai = (1 + - l)g)(l + {n- l)q') - {q, + {n - 1)%)^ (23) 
Ao = (l-g)(l-g')-(9i-go)'. (24) 
Substituting these terms into ( [T9l) . we get 

0Rs(/3, n) = n^{f3' (1 + (n - l)gP) + (fS*)' {I + {n - l)(g')^) 



4 

+ ^ (In Ai + (n - 1) In Aq) +n{l + ln(27r)} . (25) 
In the limit n — )■ 0, we have 

Ai ~ Ao + n {g(l - q') + g'(l - q) - 2go(gi - go)} • (26) 
Then, we get gnsW) = lim„_^o 0rs(/3, ^)/2n as 

9Rs{P) = y (1 - /) + (1 - (g'n +m' io! - 

1 fj^^^ ^ 9(1 - g') + g'(i - g) - 2go(gi - go) 



4 1 Ao 
+ i{l + ln(27r)}. (27) 
The saddle-point conditions yield 

/^p?''"' - T2 {^(1 - + (^1 - go) (-2go + q'iqi + go))} = 0, (28) 

'^O 

- T2 {^'(1 - + - ^o) (-2go + g(gi + go))} = 0, (29) 

^0 



- ^{gi ((1 - g)(i - g ) + (gi - go)') 



- (gi - go) (g + g' - 2gg') - 2(gi - go)^ = 0, (30) 



Mo'' - i {go ((1 - g)(i - g') + (gi - go)') - (gi - go) (g + g' - 2gg')} = o, (3i) 

-^0 

where we put fip = p/3^J^/2. 



3.1.1. Remarks for RS solutions In figure [T], we assumed that a set of spin 
configurations consisting a pure state with the weight e'^'^^^^Z also consists a pure 
state with the conjugate weight (e~^^^^"*)* /Z* . This can be accepted by considering that 
the two weights of an identical spin configuration S = S' take the same absolute value 
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|g-/3W(S)^2'| = |(e^^^'^'^))*/Z*|, which imphes that the support of each weight becomes 
identical. This indicates that the above assumption holds, since each pure state can be 
regarded support of the Boltzmann weight. 

Each pure state a has its own partition function Za, and the total partition function 
is given hy Z = Za- Using this notation and the standard description of pure 
states [7j and focusing on the P2 solution as an example, we can write the physical 
significance of the overlap q as 

1 Tr 5,e"'^^(^)5,(5) Tr 5ie-^«(^)5b(5) 



N 



(32) 



where we put Wa = Za/Z and introduce an indicator function 6a{S) which is defined 
as Sa{S) = 1 if S belongs to the pure state a and Sa{S) = otherwise. Similarly, the 
conjugate overlap q' becomes 

1 TrS,(e-'5«(^))*(5,(5) Ti Si{e-^'^(^^y 6^(8) 



N 



and the inter-overlaps qi and go are expressed as 

1 TiSiC-^^^^'^daiS) Tr5,(e-''«(^))*5,(5) 



qi = Y^ Wawl 



N 



% = Y^aw, Y — z: ^ • 

aj^b i " ^ 



(33) 

(34) 
(35) 



Combining (I32l)-( l35|) and the fact that the pure states are common for Z and Z*, we 
can derive the constraints for the overlaps as 



q =q 



qi = Qi , qo = Qo 



(36) 

We can easily find that fl28l) - fl31l) actually have the solutions satisfying these constraints. 
For the rest of this paper, we assume fl36|) from the beginning of analyses. 



3.2. IRSB 

Next, we derive the IRSB solutions of ( |T9|) . In a similar way as the RS case, we have 
two types of solutions represented in figure [2j The corresponding overlap matrices are 
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After some calculations, we get 

01RSB (/5,n,m 



n 
2p 



(1 + (m - l)q^ + {n- m)ql) 



Q 



( Qi Qi 
Qi Qi 



Qo Qo 
\ Qo Qo 



Qo Qo \ 
Qo Qo 



Qi Qi 
Qi Qi ) 



(37) 
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Figure 2. Two types of IRSB solutions. 
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+ /i; (1 + (m - l){q[Y + {n- m){q',Y) + 2|/ip| [mql + {n - m)ql) 

1 / / 71 \ ft \ 

+ - lnr/2+ 1 lnr/i + — (m-l)lnr/o + n {1 + ln(27r)} , (38) 

2 \ \m / m / 

where 

r/o = (l-gi)(l-gl), (39) 
r^i = (1 + (m - l)gi - mgo)(l + (m - l)gi - mgg) - m^(gi - go)^, (40) 
?72 = (1 + (m - l)qi + (n - m)go)(l + (m - l)gi + (n - m)go) - (mgi + (n - m)go)^ 
~ + n {go (1 + (m - l)gi - rnq'^) + g^, (1 + (m - l)gi - mgo) - 2mgo (gi - go)} • (41) 

For simplicity, let us hereafter assume the condition qo = q'o = qo = 0, which is expected 
to be satisfied due to the spin-reversal symmetry. Under this condition, we obtain 

gmsB{f3,m) = \imn^o(pmsB{(3,n,m)/2n as 

1 



9iRSBil3,m) 



Ap 



fip (1 + (m - l)qf) + (1 + (m - l){q[y) + 2m\fip\f, 



+ {Inr^i + (m - 1) Inr^o} + ^ {1 + ln(27r)} . 
4m 2 



Taking the variation with respect to gi, gi and m, we get 

= 0, 



(42) 

(43) 
(44) 

VI ^-'/oJ=0, (45) 

where the saddle-point condition with respect to q[ is omitted since it gives the complex 
conjugate of gi, the reason of which is the same as explained in section [3. 1.11 



1 / 1 + (m - l)g; 1 - q[ 
M + — i 

l/^pl^r^ - — = 0, 

^ W+/^;(g'ir+2KigT}-^ 

^ 1 / gi(l + (m-l)g^) + g^(l + (m-l)gi)-2mg,^ 



{In ?7i + (m — 1) In 770} 
+ ln ?7o 



Zeros of spherical spin glasses and chaos 



12 



4. Phase diagrams and DOZ 
4-1- p = 2 case 

It is known tliat tlie RS solution is sufficient for the p = 2 case. In this case, fl28|) and 
( 13T|) are low-degree polynomial equations and can be analytically solved. Based on the 
physical descriptions in figure [H we get the following three solutions: 

PI This solution is the usual paramagnetic solution q = q' = qi = qo = 0. The 
generating function g{l3) becomes 

= ^(/^P + /^P + ^ {1 + 1^(2^)} = - PI) J' + ^ {1 + 1^(2^)} , (46) 

and the corresponding DOZ is ppi = 0. 

P2 This solution is given hj q = q' = % = and qi > 0. Assuming g = g' = 0, we can 
easily solve (130|) as 

The inter-overlap gi should be real as explained in section I3.1.H which means that 
this solution is valid only for \f3\J > 1. Substituting this solution into (!27|) . we get 

9P2{f3) = ^ - ^ {1 + HWI'J')} + ^ {1 + ln(2vr)} . (48) 
The corresponding DOZ yields a finite value 

PP2 = ^, (49) 



which is the same value as the REM's one \i6[ [21] . 

RSSG We impose g, g' 7^ and gi = go = g in (!27j) and find that g{f3) does not depend 
on g. This is because the contribution of g is proportional to O(n^) in (IT^ and 
vanishes in the limit n 0. This means that we cannot distinguish the solutions 
RSSGl and RSSG2 in figured! and hence we just call this RSSG. The solution of 
(|28|) and ([29D is given by 

g = 1 - q' = q*. (50) 

A condition 3f?g > required in taking the spin trace of (ITTj) leads tog = l — 1//3J 
for Pi > and g = 1 -|- 1//3J for /3i < 0. The generating function then becomes 

giP) = lAI J - 1 {3 + ln(|/3|V2)} + ^ {1 + ln(27r)} . (51) 

We can easily find that the DOZ also vanishes as the PI case, Prssg = 0. 

Summarizing these results, we can derive the phase diagram and show it in the left panel 
of figure [31 The phase boundaries are obtained by equating the generating functions 
g{f3) of the adjacent phases. Note that gp2 is always larger than gpi and g-RssG except 
for on the phase boundaries, which seemingly implies that the P2 solution dominates 
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Figure 3. Complex-temperature phase diagrams of tlie p-body spherical SG for p — 2 
(left) and p = 3 (right) with J = 1. The DOZ of the P2 phase for p = 2 stays constant 
and is shaded, while that for p = 3 is not constant and is coloured depending on 
the values. The SG and PI phases have no DOZ in both the cases. On bold phase 
boundaries, the P2-P1 and P2-SG ones in the right panel, the one-dimensional DOZ 
becomes finite. The dashed line in the right panel is the extension of the boundary 
between PI and P2 phases and is not the true phase boundary. 

all the complex /3 plane, if we blindly follow the saddle-point method. Obviously, this 
is incorrect. The constraint qi = ql means the failure of the P2 solution in a region 
1/3 1 J < 1, which explains the emergence of the PI solution. On the other hand, we have 
no direct reason to explain the phase transition from P2 to RSSG. We choose the RSSG 
phase in the region /3i > 1 based on the physical appropriateness. This point will be 
clearer by considering the region n > 0, but it is beyond our current purpose in this 
paper. 

The DOZ on the boundaries should be evaluated separately by using ([6]). For 
the P1-P2 boundary, we can calculate the DOZ by putting gi = gpi, g2 = gp2 and 
b{/3i, (32) = 1 — {f3f + in IQ. Note that we cannot put b = gpi — gp2 in this case, 

since the P2 solution is always larger than or equal to the other solutions as explained 
in the above paragraph. Simple calculations show that the density vanishes. Similarly, 
we can evaluate the density on the P2-SG boundary and again it becomes zero. 

Based on the derived DOZ, we can directly calculate the specific heat C{f3) by the 
formula 

f 13"^ 
C{(3) = - J dzidz2p{zi,Z2)j^^—^. (52) 

We find that the result completely agrees with the known values [50], which validates 
our calculation. 

We also notice that the DOZ inevitably vanishes in the RSSG phases. This is 
because the absence of contributions from the inter-overlaps q, where the generating 
function g[(3) becomes analytic with respect to /3. This consideration assures that there 
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is no kind of transitions or chaos effects in the RSSG phase. 



4-2. p = 3 case 

We first derive the phase diagram on the complex-temperature plane. Here, we 
investigate only the p = 3 case, but the result is expected to be essentially common 
for p > 3. In the following discussion, we treat the 1RSB2 solution described in figure [2] 
as the correct SG solution in this case. We find that the other IRSB solution, the IRSBl 
one with gi = go = in figure [21 actually exists but it shows an unphysical behaviour. 
Hence we reject it. 

The solutions PI and P2 are essentially the same as the p = 2 case. The difference 
between gpi and g'P2 depends only on and the boundary between these phases 
becomes a circle whose radius is obtained by comparing the values of gpi and gp2 with 
substitution of the solution of (|30|) under the condition q = q' = 0. The resultant radius 
becomes \(3\ = (3p ^ 1.39884/ J. 

The SG phase of the p = 3 case is known to be described by the IRSB solution [5T] . 
We derive the transition temperature Pc from PI to SG at the real axis in our 
formulation. The equations of state f H3|) - fH5l) have some solutions even at (32 = 0. 
A reasonable solution among them is obtained under the condition qi = q'l = qi, since 
the usual overlap matrix with a real temperature is recovered by this condition. Actually 
from (l43!) - fH5l) . we can derive the same equations of state as the one under the usual 
real-parameter case by putting qi = q'l = qi and assuming (3 is real |5l], though the 
breaking parameter m in the usual case is replaced by 2m in our formulation. This is 
natural since each replica is doubled to calculate = (ZZ*)" in our formulation. 

Hence, additionally assuming the condition m = 1/2 (which corresponds to m = 1 in 
the usual formulation), we can calculate the transition temperature (3c by solving (H3ll 
and (l45l) . which leads to the known value (3c ~ 1.70633/ J [5T] . 

We can easily confirm that (/irsb with the condition m = 1 accords with gp2- This 
implies that the boundary between P2 and SG phases is obtained by solving ( H3l) -( l45l) 
under the condition m = 1, which is actually the case for the REM [21j. This can also 
be seen from that the IRSB equation of qi ( l44l) coincides with that of P2 ( l30l) at m = 1. 
Solving ( H3l) - (H5ll under the condition m = 1 involves some technical difficulties, the 



details of which are given in [Appendix A 



In the SG phase, we need to directly treat all the equations of state (I43l)-(l45l) to 
calculate giRSB- The phase boundary between the SG and PI phases is obtained by 
equating gpi and gmsB- The technical difficulties to evaluate (H3l)-(H5l) in this case are 



also summarized in Appendix A 



Summarizing the above points, we can write the phase diagram for the p = 3 case 
and give it in the right panel of figure [3l We can find some difference from the p = 2 
case in the shape of the diagram. The shape is actually related to the DOZ on the 
boundaries, which is explained below. 
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4-2.1. Density of zeros Here, we calculate the DOZ for the p = 3 case. 

The PI phase is the trivial case. We can easily confirm that the DOZ of PI phase, 
Ppi, is uniformly zero as for p = 2. 

The case of the P2 is more complicated. We can rewrite the generating function as 
gp2 = gpi + A(/3i,/32,gi(/3i,/32)), where 

A(/3i, /32, P2)) = + Pl)J'q'i + \ ln(l - qf), (53) 

and the finite contribution to pp2 only comes from A. Differentiating A with respect to 
/3i twice, we get three terms 



d'^A\ fd^A\ f d^A \ fdqA fd^A\ f dq 



where we omit a term being proportional to dA/dqi since it vanishes due to the saddle- 
point condition. Subscripts of the brackets denote the fixed variables in taking the 
partial differentiation. Evaluating each term yields 

'd^A\ 1 



d"^^ \ Po i2~p-l 



(/3i' + /32')Ar -^TIT^' (57) 

and the factor dqi/dPi can be calculated by differentiating (!30|) with respect to The 
result is 



dfi, l2gi-K|(p-2)gr(l-gf^2 



The counterpart with respect to /32 can be obtained as well. Summing up both the 
contributions, we get 

J^~v , p^(/3? + /3|)At-^y 
= 2^^^ + 2^^ 

2/02 I /o2\ 74;;2p-4 ^ 1 i ;;2 



We can easily evaluate this equation after solving ( l30l) . In figure [3l pp2 is coloured on 
the phase diagram. We can find that the density tends to decrease as P>2 grows. To see 
the quantitative behavior more precisely, we plot the density on the imaginary axis in 
figure m The figure shows that the value of the density converges to a constant as 
grows, the limiting value of which is J'^ jln as the p = 2 case, which can be understood 
from dSni) and dSH]). 

The DOZ of the SG phase can be assessed in a similar way to the P2 case, although it 
requires rather involved calculations due to the existence of four variational variables qi , 
q[, qi and m. The resultant formula is not enlightening and we here omit it. Evaluating 
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Figure 4. The DOZ on the imaginary axis (3i — for the p = 3 case with J = 1. 
The vertical dashed line denotes the transition point /32 — Pp between the PI and P2 
phases and the horizontal one represents a constant J^/27r being the limiting value of 
PP2 in the limit j32 — > oo. 

the DOZ through the formula, we find that there are no zeros in the SG phas^, which 
is the same as the REM. This imphes the absence of any kind of phase transitions in 
the SG phase of this system. 

Next, we refer to the DOZ on the phase boundaries. First derivatives of g{[i) of 
adjacent phases are crucial as already noted in (El). 

On the P1-P2 boundary, choosing gi = gpi, g2 = gp2 and b = gpi — gp2 in ([6]), we 
find the one-dimensional DOZ ppi_p2 as 

pPi-P2(/3i, P2) = \j'q? {P! + Pi) S (A(/3i, /32, qi)) , (60) 

where A is given in (1531) . Clearly, this yields nonzero value, which is natural since the 
transition between the PI and P2 phases is of first order in the p = 3 case. Note that 
the derivative of qi, which involves the factor dgp2/dqi, does not appear in this formula, 
since dgp2/dqi vanishes due to the saddle-point condition. 

We can also confirm that the DOZ takes finite values on the Pl-SG boundary by 
simple calculations. Since there are no zeros in both the phases, the Yang-Lee theorem, 
which proves no phase transition in a region without the DOZ, requires that the one- 
dimensional DOZ on the boundary cannot vanish. This is in contrast to the p = 2 case 
where the P2 phase intercepts the PI and SG phases. 

Meanwhile, on the P2-SG boundary, the one-dimensional DOZ does not appear, 
since the condition m = 1 at the boundary makes the first derivatives of gp2 and giRSB 
identical. This observation implies that there is no one-dimensional density on the P2- 
SG boundary in general, which is expected to be applicable to other situations and other 
SG systems. 

Before closing this subsection, we mention the physical consequence of the absence 
of the zeros in the SG phase in the present situation. From the point of view of the 

§ We numerically evaluate the DOZ at several regions in the SG phase, such as around the P2-SG 
boundary and on a line /3i = /3c, and confirm that the values are smaller than 10~^^ at most. This is 
in a margin of numerical errors of our calculation and we conclude that the DOZ is zero. 
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chaos effect, this is quite natural since the present model does not show the temperature 
chaos [30] . Generally, the chaos effect is connected to the change of the dominating pure 
states when we vary the corresponding physical parameter. For the spherical SG model, 
it is shown that the pure state essentially does not change in temperature below the 
SG transition point. This means the absence of the temperature chaos, and of a certain 
transitions mentioned in section [1] To provide a clear connection between the chaos 
effect and the DOZ, we need to investigate models exhibiting the chaos effect. One of 
the simplest choice of such models is the spherical SG model with multiple many-body 
interactions [30], which is examined in the next subsection. 

4-3. Multiple interaction case 

The Hamiltonian of the {p + r)-body interacting spherical model is given by 

H = — ^ ^ Ji^...i^Si-^ ■ ■ ■ Si^ — t ^ ^ Ki-^...i^Si-^ ■ ■ ■ Si^. (61) 

ii<---<ip n<-<ir 

As ([9]), the coupling constants Ji^-i^ and Ki^...i^ are Gaussian variables with the variances 
.Pp\/2NP~^ and JM/2N'-^, respectively. 

The calculation of (f){P,n) of this system is the same as given in section [31 The 
resultant expression of (f){/3,n) is obtained by just replacing g^j^ with (<?„;, + e^g^^) in ( fT9l) . 
and {qaf,y and q^f^ are replaced in a similar manner. The analysis to obtain the phase 
diagram and the DOZ is also the same and we omit it here. We just give the result 
below. 

As an example, we show the case {p,r) = (3,4) and e = 0.2, where the IRSB 
solution is correct [30] . The shape of the phase diagram is almost the same as that of 
the usual p = 3 spherical SG model, and the corresponding critical temperatures are 
Pc = 1.68462/ J and Pp = 1.37877/ J. A remarkable property in this case is that the DOZ 
in the SG phase is finite. In figure [5], we give the DOZ coloured in a logarithmic scale on 
the phase diagram. We can find that the DOZ values of the SG phase are rather small 
in comparison with the P2 one: The typical values are p = O(10~^) and p < O(10~^) in 
the P2 and SG phases, respectively. This means that the DOZ discontinuously changes 
at the P2-SG boundary as for the single p case. 

Obviously, the DOZ on the real axis is more important, and we plot it in the left 
panel of figure [6l The figure clearly shows that the DOZ on the real axis is finite. This 
is quite contrast to the case of the single p-body interaction, which strongly suggests 
that the DOZ on and around the real axis in the SG phase is closely related to the 
temperature chaos. We further discuss about this point in the next subsection. For 
quantitative comparison, we also give the DOZ on a vertical line at the critical point 
/3 = /9c in the right panel of figure [6l From these figures, we can see that the density tends 
to decrease/increase as the real/ imaginary part of f3 grows. Although it is interesting 
to see the limiting value of p in the limit /32 — oo in the right panel of figure |6l we 
could not find appropriate solutions for large (32 due to some uncontrollable technical 
difficulties. 
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Figure 5. Complex-temperature phase diagrams of the {p + r)-body spherical SG for 
p = 3 and r = 4 with J = 1 and e = 0.2. The DOZ takes finite values on the P2 and 
SG phases and is coloured in a logarithmic scale. The one-dimensional DOZ becomes 
finite on bold phase boundaries as the single p = 3 case. 
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Figure 6. The DOZ on the real axis (left) and on a vertical line at /3i = f3c (right) in 
the SG phase of the case {p, r) = (3, 4) and e — 0.2. 



Discussion 

So far, we have investigated the zeros of the spherical SGs in several situations. For 
the p = 2 case where the RS solution is correct in all the region, the zeros distribute 
two-dimensionally in the P2 phase, though the SG phase has no zeros. The zeros are 
also absent in the IRSB-SG phase for the p = 3 case. These observations indicate that 
the zeros are not directly related to the RSB. On the other hand, in the SG phase of 
the (p -|- r)-body interacting spherical model, the DOZ takes nonzero values. Since this 
system shows the temperature chaos being absent in the single p-body interaction case, 
we can naturally speculate that the chaos effect is closely related to the DOZ distributed 
in the SG phase, especially the values of DOZ on the real axis is quite important. 
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The above statement becomes clearer by referring the DOZ of the REMs |21j . 
For the generahzed REM (GREM) in the continuum hmit of the hierarchy, the one- 
dimensional DOZ on the phase boundaries are accumulated to become a two-dimensional 
distribution. That two-dimensional DOZ looks similar to the (p-l-r)-body spherical case, 
but a crucial difference is that the DOZ of the GREM is zero on the real axis of (3. This 
can be interpreted as follows. In the GREM, by definition, the consecutive transitions 
in temperature are characterized by a sequential freezing of a part of the total spins |54j . 
On one hand, this means that the nature of the transitions essentially becomes of second 
order, which explains the absence of zeros directly on the real axis. On the other hand, 
the freezing character of the transitions implies the presence of correlations among the 
equilibrium spin configurations at different temperatures, which leads to the absence 
of the temperature chaos as shown in [26]. Combining this observation and the fact 
that the chaos effect is connected to discontinuous changes of the dominant part of pure 
states [30] , we can reasonably conclude that the chaos effect generally has a first-order- 
transition-like nature and is signaled by the zeros on the real axis. 

We here refer to the Sherrington-Kirkpatrick (SK) model [55] known to show the 
temperature chaos |29]. In this model, although it is known that there emerge many 
pure states in the SG phase, properties of the pure states are quite unclear in comparison 
with those of the spherical models. This may be regarded as a consequence of a local 
instability of the phase space: The so-called de Almeida-Thouless instability [56] . This 
instability directly causes the SG transition of the SK model, and accordingly the nature 
of the transition becomes continuous, which may question the above description of the 
chaos effect. However, considering Bethe SGs expected to show the essentially same 
behavior as the SK model, we can find that the zeros distribute in the SG phase even 
on the real axis of temperature pQ] . This fact implies that the discontinuous change of 
the dominant pure states also occurs in the SG phase of the SK model, which supports 
our above description. 

We should also notice that the zeros on the real axis in the SG phase do not mean the 
singularities of the free energy. Actually, the free energy of this model for real /3 does not 
show any singularities in the SG phase. This is because the DOZ are two-dimensionally 
distributing in the SG phase. An interpretation of this fact can be naturally obtained by 
considering the analogy with the electrostatics explained in section 12.11 We can see that 
the electrostatic potential can be analytic in a region where the point charges distribute 
two-dimensionally, which also means the absence of the singularities of the free energy 
in a region where the zeros distribute two-dimensionally. 

Conversely, partition-function zeros can detect extraordinary behaviour even not 
appearing as singularities of the free energy. We stress that our result in this paper 
becomes the first evidence and demonstration of that fact. 
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5. Conclusion 

In this paper, we investigated the zeros in the complex temperature plane of the many- 
body interacting spherical SGs, in the single p-body and multiple (p+r)-body interacting 
cases. Our formulation utilizes the replica method and generalizes the Parisi scheme 
to be applicable in the complex-parameter case. The relations between the pure-state 
structure and the overlap matrices were also considered. Based on the formulation, we 
derived the phase diagrams in the complex-temperature plane and calculated the DOZ 
on each phase and boundary. By changing the parameters p and r, we could easily 
investigate several physically-different situations and examined the possible relations 
among the DOZ, the RSB and the chaos effect. The notable significances of the result 
are as follows: 

• The RSSG phase cannot have the finite DOZ. 

• The RSB is not necessarily connected to the two-dimensionally distributed zeros. 

• The two-dimensionally distributed zeros around and on the real axis are closely 
related to the chaos effect, and do not necessarily lead to singularities of the free 
energy. 

The formulation presented in this paper has some possible applications. One of the 
most simple applications is to investigate the SG systems with the FRSB such as the SK 
model. Although we constructed our solution in the IRSB level, it is possible to extend 
the solution to the FRSB. The DOZ in the FRSB phase of the Bethe SG was studied 
in [20], but the RS ansatz was used to calculate the DOZ, which can involve potential 
errors in the estimation. Besides, in that result we cannot see the discrimination between 
the P2 and FRSB-SG phases, which is clearly different from our present result. This 
difference may be due to the RS ansatz [20], or due to the difference between the IRSB 
and FRSB. It is also a question how the de Almeida- Thouless instability [56] relates to 
the DOZ. These questions motivate us to investigate the FRSB-SG phase by our current 
formulation, which will be our future work. 

In the presented paper, we extracted some physical properties of the zeros only 
on and around the real axis, but the DOZ in the whole complex plane potentially 
has more significance. The P2 phase has a particular interest since the DOZ on and 
around the imaginary axis is possibly related to some dynamical properties [6], l57] . 
Hence, it will be interesting to seek physical contents of the DOZ in the whole complex 
plane. Other possible applications, such as to finite-dimensional SGs [32] and structural 
glasses [58], are also important directions. We hope that the presented formulation and 
result inspire those researches, and lead to revealing origins of extraordinary behaviour 
of glassy systems. 
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Appendix A. Solving the equations of state with complex parameters 

Here, we describe how to solve the complex equations of state and to obtain the phase 
boundaries in figure [31 

We start from the simpler case, i.e. the boundary between the P2 and IRSB-SG 
phases. As explained in section 14. 2[ the condition m = 1 is essential to determine 
the P2-1RSB boundary. Under this condition, the generating functions gp2 and Qirsb 
automatically becomes identical. This means that Qirsb becomes independent from 
qi and q[, and the same is true for the equation of state of qi, (jH]). On the other 
hand, the variational condition with respect to m, ( l45l) . which should be satisfied at 
the phase boundary, still depends on qi and q[. Hence, to obtain the phase boundary 
we need to calculate qi through (H3|l . though it is not needed to evaluate (^irsb at 
m = 1. Summarizing these observations, we adopt the following procedures to obtain 
the boundary between the P2 and IRSB-SG phases: 

(i) Fix a value of /i = p/3^J^/2. 

(ii) Calculate gi through (jH]) with m = 1. This can be easily performed by usual 
methods such as iteration, bisection method and Newton's method. For p = 3, 
even the analytic solution can be obtained. 

(iii) Calculate qi by solving f H3|) with substitution of the above qi under m = 1 (note 
that q'l = ql )■ Since qi is complex, the iteration and bisection methods do not 
work well. We employ Newton's method with an appropriately-chosen initial value 
of qi. Empirically, we find that the initial value should have a small imaginary part 
and a real part slightly smaller than unity. 

(iv) Evaluate the left-hand side of ( H5l) with m = 1 by using the obtained values of qi 
and qi. If the value is sufficiently small, the given /3 = ^J^iifpl J gives the desired 
phase boundary. Otherwise, restart from (i) with a new value of /x. To efficiently 
search the boundary, we actually fix 53/3 and gradually change 3fJ/9 to find a value 
of /3 at which ( H5|) is satisfied. 

These procedures actually work well and the P2-1RSB boundary is obtained 
straightforwardly. 

Next, we consider the boundary between the PI and IRSB-SG phases. For this, 
we need to treat all the equations of state ( l43|) -(H5l) with mutual dependence among the 
order parameters, unlike the P2-1RSB case. To actually solve ( 143|) - (H5|) . we focus on the 
fact that q\ and the left-hand side of ( H5ll are real. This property enables us to use the 
bisection method to evaluate those two equations. The resultant procedures we accept 
are as follows: 



(i) Fix a value of /i = p/3^J^/2. 
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(ii) Calculate qi, qi and m. Call the bisection subroutine with respect to m by fl45l) . 

(a) Fix three appropriate values of m, > m' and = {m^ + m')/2. 

(b) Calculate qi and qi for given three values of m. For this, call the bisection 
subroutine with respect to qi by (jH]). 

1. Fix three appropriate values of qi, q^ > q{ and = (gj* + q{)/2. 

2. Calculate qi for given m and three values of qi by solving (H3|) and obtain 

g™5 Q'l)- Newton's method with an appropriate initial value of qi is 
again useful. 

3. Compare the left-hand side of (jHj) with substitutions of {qi,qi), {q^^AT) 
and {q\,(fi). Replace q^ or q{ with depending on the compared values 
of dH. 

4. Repeat 1-3 until the value of q converges. Return the resultant (gi,gi). 

(c) Compare the left hand side of (jH]) with substitutions of {m^,qi,qi), 
{m^,q^,q^) and {m\q\,q\). Replace or m} with depending on the 
compared values of ( H5i) . 

(d) Repeat (a)-(c) until the value of m converges. Return the resultant (m, gi, gi). 

(iii) Compare the values of g^i and Qiksb by using the obtained (m,gi,gi). If the 
difference is sufficiently small, the given gives the phase boundary. Otherwise, 
restart from (i) with a new value of /i. 

To actually conduct these procedures, some difficulties are involved in choosing the 
appropriate values of gi and m for the bisection subroutines. If the chosen values are 
inappropriate, the converged values become unphysical. To resolve this point, we start 
from the critical point (3 = (3c and gradually change the value of (3 step by step, with 
using the values of order parameters in the previous step as the initial values of gi and 
m in the current step. Although this prescription works well for assessing the phase 
boundaries, the trouble becomes more serious when we evaluate the DOZ in a region far 
from the real axis in the SG phase, due to bad behaviour of ( I43l) - (l45l) in that region. A 
more effective routine to solve the complex equations of state (H3l) -( H5l) will quite benefit 
to assess the DOZ, but in the presented results we did not pursue this point and just 
tuned the parameters until physically-plausible results are obtained. 
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